##### 加载包、函数
# 加载相关包
library(openxlsx)
library(reshape2)
library(geepack)
library(gtsummary)
library(officer)
library(flextable)
library(mice)

# 加载自定义函数
source("./func/gee_run.r")

# 用一个元素为 gee_fit 的二维 list() 绘制所需的表格
summary_out = function (fit_table) {
  # summary_out
  ret = list()
  for (i in 1:length(fit_table[[1]])) {
    tb1 = list()
    for (j in 1:length(fit_table)) {
      tb1[[j]] = tbl_regression(fit_table[[j]][[i]], exponentiate = T, include = c(1))
      print(i)
    }
    ret[[i]] = tbl_merge(tb1, tab_spanner = varsy)
    
  }
  return (ret)
}
#####


##### 包配置
# 设置 docx 的样式
sect_properties = prop_section(
  type = "nextPage",
  page_size = page_size(),
  page_margins = page_mar()
)
#####


###### 读取数据，数据预处理
data = read.xlsx("./input/data1105.xlsx", 1)
varsx = colnames(data)[c(14:25, 27, 65)]     # 自变量
varsy = colnames(data)[c(9:13)]              # 因变量
varsc = colnames(data)[c(31:32, 34:39, 66)]  # 协变量
# 保留所需的列
data_test = data %>%
  select('id', 'month', varsy, varsx, varsc)

# 筛选：8～24 月龄中填写次数 >=2 的个体的所有观测，然后去掉 month == 6 的观测
extract_id = data_test %>%
  dcast(id ~ month) %>%                          # 重铸后的变量列表为 c(`6`, 8`, `12`, `18`, `24`)
  mutate(count = `8` + `12` + `18` + `24`) %>%   # 添加一列 count，计算每个个体 8～24 月龄填写的次数
  .$id                                           # 得到符合条件的个体 id
data_anal = data_test %>%
  filter((id %in% extract_id) & month != 6)      # 取得这些个体的所有观测（但不要 6 月份的观测）

# 确保变量类型一致：因变量为数值，自变量、协变量为因子
data_anal$id = as.character(data_anal$id)
data_anal[varsy] = lapply(data_anal[varsy], as.numeric)  # gee 要求因变量是数值
data_anal[varsx] = lapply(data_anal[varsx], as.factor)
data_anal[varsc] = lapply(data_anal[varsc], as.factor)

# 删除自变量缺失的行
data_anal = data_anal[complete.cases(data_anal[, varsx]), ]
data_anal[varsx] = lapply(data_anal[varsx], as.factor)      # 删除后，有些因子的计数可能为 0，需要重置 factor
######


##### mice
data_tmp = unique(data_anal[, c("id", varsc)])  # 只填补 varsc 的列
miss = md.pattern(data_tmp)                     # 查看缺失数据的模式
imp = mice(data_tmp, m = 5, seed = 1)           # m=5 指定 mice 填补 5 组数据
varsc_mice = complete(imp, action = 1)           # action=1 指示选择第一组数据
data_mice = merge(data_anal[, c("id", "month", varsy, varsx)], varsc_mice)   # 替换原来的 varsc 区域数据
write.xlsx(data_mice, './output/1105_mice_data.xlsx', overwrite = T)         # 保存以便后续使用

# 查看 mice 前后的 varsc 分布情况 - tbl_summary
varsc_before_mice = data_anal %>%
  select('id', varsc) %>%
  unique() %>%
  select(varsc) %>%
  tbl_summary()
varsc_after_mice = data_mice %>%
  select('id', varsc) %>%
  unique() %>%
  select(varsc) %>%
  tbl_summary()

save_as_docx(
  values = lapply(list(varsc_before_mice, varsc_after_mice), as_flex_table),
  path = "./output/1105_varsc_sumary_compare.docx",
  pr_section = sect_properties
)
# 可以用以下代码直接查看对比图
# tbl_merge(tbls = list(varsc_before_mice, varsc_after_mice), tab_spanner = c('before', 'after'))
#####


##### gee，并使用 tbl_regression 定制表格输出
data_mice = read.xlsx("./output/1105_mice_data.xlsx", 1)
data_mice[varsy] = lapply(data_mice[varsy], as.numeric)
data_mice[varsx] = lapply(data_mice[varsx], as.factor)
data_mice[varsc] = lapply(data_mice[varsc], as.factor)

fit_list = mice_data %>%
  gee_run(
    varsx = varsx,
    varsy = varsy,
    #varsc = varsc,
    file = "./output/1105_mice_data_gee.xlsx"
  ) %>%
  summary_out()


# 保存到 docx
save_as_docx(
  values = lapply(fit_list, as_flex_table),
  path = "./output/1105_mice_data_gee_tbl_regression.docx",
  pr_section = sect_properties
)
#####
